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Modulated Oscillations in Three Dimensions 

Jonathan M. Lilly, Member, IEEE 



Abstract — The analysis of the fully three-dimensional and time- 
varying polarization characteristics of a modulated trivariate, 
or three-component, oscillation is addressed. The use of the 
analytic operator enables the instantaneous three-dimensional 
polarization state of any square-integrable trivariate signal to be 
uniquely defined. Straightforward expressions are given which 
permit the ellipse parameters to be recovered from data. The 
notions of instantaneous frequency and instantaneous bandwidth, 
generalized to the trivariate case, are related to variations in the 
ellipse properties. Rates of change of the ellipse parameters are 
found to be intimately linked to the first few moments of the 
signal's spectrum, averaged over the three signal components. In 
particular, the trivariate instantaneous bandwidth — a measure 
of the instantaneous departure of the signal from a single pure 
sinusoidal oscillation — is found to contain five contributions: 
three essentially two-dimensional effects due to the motion of the 
ellipse within a fixed plane, and two effects due to the motion of 
the plane containing the ellipse. The resulting analysis method 
is an informative means of describing nonstationary trivariate 
signals, as is illustrated with an application to a seismic record. 

Index Terms — Instantaneous frequency, instantaneous band- 
width, nonstationary signal analysis, trivariate signal, three- 
component signal, polarization. 

I. Introduction 

MODULATED trivariate or three-component oscillations 
are important for their physical significance. A wide 
variety of wavelike phenomena are aptly described as modu- 
lated trivariate oscillations, including seismic waves, internal 
waves in the ocean and atmosphere, and oscillations of the 
electric field vector in electromagnetic radiation. Real-world 
waves often appear as isolated packets, as evolving nonlinear 
wave trains, or as sudden events whose properties change 
with time — all situations involving nonstationarity. To date 
interest in trivariate signals has primarily been motivated by 
seismic applications. In oceanography, measurements of the 
three-dimensional velocity field have traditionally been rare, 
but recent improvements in both measuring and modeling the 
three-dimensional oceanic wave field, as in (TJ, ||2) and (3), 
10) for example, make the analysis of trivariate oscillations 
increasingly relevant to this field as well. Therefore a suitable 
analysis method for nonstationary trivariate signals would find 
broad applicability across a variety of disciplines. 

Complex-valued three-vectors, introduced by Gibbs in 1884 
0, have long been used to describe the polarization state 
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of an oscillation in three dimension^. Previous signal anal- 
ysis works have considered the three-dimensional but time- 
invariant polarization of trivariate signals IfTI- lHTI . as well 
as the time-varying two-dimensional polarization of bivariate 
signals, potentially along a set of three orthogonal planes [ 12 1- 
lfl6| . The purpose of this paper is to enable the analysis of the 
fully three-dimensional instantaneous polarization state of a 
modulated trivariate oscillation, and to relate the variability of 
the polarization state to the moments of the signal's Fourier 
spectrum. This is a natural but non-trivial extension of recent 
work on modulated bivariate oscillations by |[T6l . 

One approach to the analysis of trivariate oscillations is 
in terms of a frequency-dependent polarization, with key 
contributions found in a series of works by Samson Q, 
ID, 03-1111. Other authors, e.g. JlO], OH, have similarly 
modeled trivariate signals as oscillatory motions with three- 
dimensional but time-invariant polarizations, possibly in the 
presence of background noise. Estimation of the frequency- 
dependent polarization state based on the multitaper spectral 
analysis method of Thomson Ell was accomplished by Q. 
There the averaging necessary to estimate the spectral matrix 
was accomplished with an average over taper "eigenspectra" 
rather than an explict frequency-domain smoothing. 

The extension to a time- and frequency-varying three- 
dimensional polarization was pursued by [23|-[27|, by em- 
ploying multiple continuous wavelets, rather than multiple 
global data tapers. However, a limitation of this approach 
is that the polarization is a function of the multi-component 
wavelet transform and not an intrinsic property of the signal; 
thus the definition of polarization is basis-dependent. The 
time/frequency averaging implied by the use of multiple 
wavelets may introduce unwanted bias into the estimate, but 
the extent of this bias is impossible to quantify because the 
polarization is not independently defined. Furthermore, re- 
liance on the wavelet basis to define time-varying polarization 
sidesteps the question as to what kind of object the signal of 
interest is, if it is in fact not a sinusoid. 

A more compelling, and ultimately more powerful, approach 
is to begin with a model of the signal itself. In the univariate 
case, the notion of a modulated oscillation is made precise 
through the use of the analytic signal |28|-[32|. This construc- 
tion permits a unique time-varying amplitude and phase pair 
to be associated with any square-integrable real-valued signal, 
see e.g. ||32l and references therein. In terms of the analytic 
signal, intuitive and informative time-varying functions may 
then be found — the instantaneous frequency |28]-|32] and 
instantaneous bandwidth If33l - ll3~5l — that formally provide the 
contributions, at each moment in time, to the first-order and 

'Gibbs referred to complex- valued three- vectors by the now- archaic term 
"bivectors", not be confused with the bivectors of geometric algebra (6). 
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second-order global moments of the signal's Fourier spectrum. 
In this way time-dependent amplitude and frequency can 
seen as properties of the signal. In noisy or contaminated 
environments, time-frequency localization methods such as 
wavelet ridge analysis [ 36 1— [ 38 1 can then be employed to yield 
superior estimates of these well-defined signal properties. 

The instantaneous description of modulated oscillations 
using the analytic signal, and the associated instantaneous 
moments, has been extended to the bivariate case by several 
authors lfT2l - |[T6l . The use of a pair of analytic signals permits 
the description of a bivariate signal as an ellipse with time- 
varying properties, as appears to have first been done by Rene 
et al. [12 1 following an application of the analytic signal to 
univariate seismic signals by |39l . It is a testament to the 
broad relevance of these ideas that there exist two distinct lines 
of development: one in the geophysics community [12], [13|, 
[39 1, and another originating in the oceanographic community 
[14 1, [16] based on an earlier body of work on the stationary 
bivariate case B51-ll44l. 

This paper extends the analytic signal approach to in- 
vestigating instantaneous signal properties to the trivariate 
case. The structure of the paper is as follows. The unique 
representation of a modulated trivariate oscillation in terms of 
a trio of analytic signals is found in Section GJJ This enables 
any real-valued trivariate signal to be uniquely described as 
tracing out an ellipse, the amplitude, eccentricity, and three- 
dimensional orientation of which all evolve in time. Simple 
expressions are derived which give the time-varying ellipse 
properties directly in terms of the trivariate analytic signal. 
In Section Hill the trivariate generalizations of instantaneous 
frequency and bandwidth are found and are expressed in terms 
of rates of change of the ellipse geometry. It is shown that five 
distinct types of evolution of the ellipse geometry can give 
identical spectra. Application to a seismic signal is presented 
in Section [TV] followed by a concluding discussion. 

All numerical code related to this paper is made freel 
available to the community as a package of Matlab routine; 
as discussed in Appendix [Al 

II. Modulated Trivariate Oscillations 

This section develops a representation of a modulated 
trivariate oscillation as the trajectory of a particle orbiting a 
time-varying ellipse in three dimensions. Unique specifications 
of the ellipse parameters are found in terms of the analytic 
parts of any three real-valued signals. 

A. Cartesian Representation 

A set of three real-valued amplitude- and frequency modu- 
lated signals may be represented as the trivariate vector 





x(t) 




a x (t) cos (j) x (t) 




x(t) = 


y(t) 




a y (t) cos 4> y {t) 


(1) 




z{t) 




_ a z (t) cos 4> z {t) 





which is herein assumed to be zero-mean and square- 
integrable. The representation ([T} is non-unique in that more 

2 This package, called Jlab, is available at http://www.jmlilly.net 



than one amplitude/phase pair can be associated with each 
real-valued signal, see e.g. [32|. However, a unique specifi- 
cation of the amplitudes a x (t), a y (t), and a z (t) and phases 
<p x (t), 4> y {t), and <fi z (t) may be found from combining x(t) 
with its Hilbert transform 



i r 

71 J-, 



x O) 

t — u 



du 



(2) 



where "f" is the Cauchy principal value integral. The six 
quantities appearing on the right-hand side of (fl~|i are taken 
to be this unique set of amplitudes and phases, called the 
canonical set, which is found as follows. 

Pairing the real-valued signal vector with i = times 
its own Hilbert transform defines the analytic signal vector 



x + (t) = 2Ax(t) = x(t) + iHx(t) 



(3) 



where A is called the analytic operator JSJJ, 11321 . The Fourier 
transform of x+(i) is given by 



e _ia, *X + (t) dt = 2U(lo)X(lu) 



(4) 



where X(w) is the Fourier transform of x(t) and U(uj) is the 
Heaviside unit step function; this follows from the frequency- 
domain form of the analytic operator. The amplitudes and 
phases of the components of the analytic signal vector 



x+(t) 



x+(t) 






y+(t) 






z+(t) 







(5) 



define the canonical set of amplitudes and phases associated 
with x(t), with a x (t) = \x+(t)\ and tfi x (t) = a.rg{x+(t)} and 
so forth; here "arg" denotes the complex argument. The real- 
valued signal vector is then recovered by x(t) — 5ft{x + (i)}, 
where "5ft" denotes the real part. 

That there is a strong physical motivation in representing 
a univariate modulated oscillation via the canonical amplitude 
and phase is now well known, see e.g. [30|-[32|, l35l . Among 
the desirable features of the canonical amplitude and phase is 
an intimate connection between these time-varying quantities 
and the Fourier-domain moments of the signal, which are 
made use of in Section [En] However, the analytic signal vector 
describes the three signal components in isolation from each 
other, whereas the fact that we have grouped these time series 
together implies that there is a reason to believe they are 
somehow related. 



B. Ellipse Representation 

Rather than consider x(t) as a set of three disparate signals, 
it is more fruitful to introduce a representation which reflects 
possible joint structure. A set of three sinusoidal oscillations 
along the coordinate axes, each having the same period but 
with arbitrary amplitudes and phase offsets, traces out an 
ellipse in three dimensions. This suggests that a useful rep- 
resentation for a modulated oscillation will be in terms of an 
ellipse having properties that evolve with time. 
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An alternate form for the analytic signal vector, the modu- 
lated ellipse representation, is therefore proposed as 



Ellipse Schematic in 3D 



x + (t)= e ^ t )j 3 (a(i))J 1 (/3(i))J 3 (^)) 



a(t) 
-ib(t) 




(6) 



where we have introduced the rotation matrices 



Ji(0) = 



J 3 (a) 




cos( 
sin t 



cos a 
sin a 






- sin ( 
cos 9 

- sin a 
cos a 
1 



(7) 



(8) 



about the x and z axes respectively. In ©, the real-valued 
signal x(t) = 5R{x + (i)} is described as the trajectory traced 
out in three dimensions by a particle orbiting an ellipse with 
time-varying amplitude, eccentricity, and orientation. It will be 
shown shortly that to each value of the analytic signal x + (t), 
one may assign a unique set of the ellipse parameters. 

A sketch of an ellipse in three dimensions with all angles 
marked is shown in Fig. Q] The interpretation of © is as 
follows. An ellipse with semi-major axis a(t) and semi-minor 
axis b(t), with a(t) > b(t) > 0, originally lies in the x— y plane 
with the major axis along the x-axis of the coordinate system. 
The ellipse is then transformed by (i) rotating the ellipse by 
the precession angle 9(t) about the z-axis; (ii) tilting the plane 
containing the ellipse about the x-axis by the zenith angle f3(t); 
and finally (iii) rotating the normal to the plane containing 
the ellipse by the azimuth angle a(t) about the z-axis. The 
position of a hypothetical particle along the periphery of the 
ellipse is specified by (j)(t), called the phase angle. All angles 
are defined over (— 7r,7r], except for f3 which is limited to 
[0, 7r], for reasons to be discussed shortly. Note that the rotation 
in three dimensions has been represented in the so-called z- 
x-z form, as this proves convenient for the subsequent matrix 
multiplications. 

One may replace the semi-major and semi-minor axes a(t) 
and b(t) with two new quantities 



K(t) 

A(t) 



a 2 (t) + b 2 {t) J_ 
2 ~ V2 



x+(t)|| 



a 2 {t)-b 2 {t) 
~ a 2 (t)+b 2 (t) 



(9) 
(10) 



the former being the root-mean-square ellipse amplitude, and 
the latter a measure of the ellipse shape; note ||z|| = V z H z 
is the norm of some complex-valued vector z, with "H" 
indicating the conjugate transpose. The quantity X(t), like the 
eccentricity y/T— b 2 (t)/a 2 (t), varies between zero for circu- 
lar motion and unity for linear motion, and thus may be termed 
the ellipse linearityU Note that yjl - \ 2 (t) n 2 (t) = a(t)b(t) 
gives the squared geometric mean radius of the ellipse, which 
could also be interpreted as a measure of the circular power, 
while X(t)K 2 (t) = [a 2 (t) - b 2 (t)] /2 could be interpreted as 
the linear power. 

3 Note ['163 defines the linearity as a signed quantity, whereas here A(t) > 0. 
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Fig. 1. Schematic of a modulated trivariate oscillation represented as 
an ellipse. A particle is shown orbiting an ellipse with constant geometry, 
characterized by semi-major axis a = 3, semi-minor axis b = 2, precession 
angle 8 = tt/3, zenith angle j3 = 7r/4, and azimuth angle a = tt/6. The 
phase increases from an initial value at 4> = 57r/6, tracing out the heavy black 
curve through one full cycle, during which time all other ellipse parameters 
are constant. The plane of the ellipse is indicated by the dotted lines, with the 
original x- and ?/-axes marked by thin dashed lines. A heavy dashed black 
line, with an open circle at its end, shows position of the "particle" at the 
initial time, while a heavy gray dashed line marks the ellipse semi-major axis. 
A heavy solid black line, with a filled circle at its end, is the normal vector to 
the plane of ellipse; the projection of this vector onto the x—y plane is shown 
with a heavy solid gray line at the top of the figure. 



The time variation of the signal is described by six rates 
of change, introduced here for future reference. The rate 
of amplitude modulation is K r (t), while A'(t) is the rate 
of deformation of the ellipse; here the primes indicate time 
derivatives. The remaining four rates of change can be said 
to be frequencies, in a sense, since they correspond to rates 
of change of angles. The orbital frequency U)$(t) = <fi'(t) 
gives the rate at which the particle orbits the ellipse. The 
orientation of the ellipse within the plane changes at the rate 
ujg(t) = 0'(t), which could be termed internal precession. This 
is distinguished from the azimuthal motion of the normal to 
the plane containing the ellipse u) a (t) = a'(t), or what we 
may call the external precession. To name the final quantity 
we may borrow a term from the description of gyroscopic 
motion and refer to uip(t) = /3'(t) as the rate of nutation; 
as this literally means "nodding", it seems to appropriately 
describe the inward or outward motion of the normal to the 
plane containing the ellipse from the vertical. 

C. Comments on Angles 

In defining the angles of the ellipse representation, we 
constrain < /3(t) < it, while the other three angles vary 
over (— 7T, 7r]. These choices deserve further comment. If the 
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ellipse geometry is constant, i.e. only <j>(t) varies in time, 
then the signal will repeatedly trace out the same ellipse in 
space, and so one should let cj>(t) vary between — tt and tt to 
accommodate such motion. Note that in (0 the substitutions 
6 n- 9 + tt and <\> M> <f> + tt both have the same effect, which 
is to change the sign of x + (t). Since <fi(t) varies between — tt 
and tt, it might appear that 8 should be limited to a range 
of tt in order to prevent this ambiguity; however, in practice 
0(t) tends to evolve continuously in situations for which the 
modulated ellipse representation is suitable, and this continuity 
means there is no ambiguity in defining 8(t) to within a factor 
of 2tt from moment to moment. 

Clearly a(t), which gives the azimuth angle of the normal 
to the plane containing the ellipse, must vary over a range 
of 2tt in order to allow for all orientations of the plane. 
The orientation of the plane can then be completely specified 
with f3(t) limited between zero and 7r/2; however, so that the 
projection of the motion onto x-y plane may be in either a 
clockwise or counterclockwise direction, f3(t) is allowed to 
vary from zero to tt. Counterclockwise motions on the x-y 
plane correspond to f3(t) < tt/2, and clockwise motions to 
(3 > tt/2. Note that this differs from the convention of [16|, 
who in their study of bivariate modulated oscillations let b(t) 
change sign to reflect the different directions of circulation 
around the ellipse. 

The Hilbert transform of x(i) decrements the phases of 
all Fourier components by ninety degrees, turning cosinusoids 
into sinusoids and sinusoids into negative cosinusoids. Thus 

Hx + (t) = -ix + (t) (11) 

by definition of the analytic signal. In the context of the 
modulated ellipse representation ©, the Hilbert transform has 
a simple geometric interpretation: the orbital phase <fi(t) of 
x + (i) is decremented by tt/2 with all other ellipse parameters 
unchanged. Thus the signal vector and its Hilbert transform 
together can be used to represent a particle moving through 
an ellipse with fixed geometry, with <fi(t) behaving as if it 
were a rapidly changing variable. This is analogous to the 
univariate case, in which the Hilbert transform of an analytic 
signal x+(t) = a x (t)e 1 ^^ shifts the phase (j) x (t) by tt/2 
with the amplitude a x (t) unchanged. 

D. Examples 

To better visualize the types of signals associated with 
the three-dimensional ellipse representation, five examples are 
presented in Fig.|2^-e; the last panel, Fig. [2J 7 , is not used until a 
subsequent section. In each of the five examples, exactly one of 
the five rates of change describing the ellipse geometry — K,'(t), 
X'(t), LUg(t), uj a (t), and uip(t) — is nonzero. As the orbital 
phase tfi(t) varies in time along with one of the five geometry 
parameters, a curve is traced out in three dimensions. The 
projection of this motion onto the x-y plane is also shown. 
The shading of the curve represents time, with the curve being 
black at the initial time and fading to light gray as time 
progresses. 

The first three examples, Fig.|2^-c, all involve the motion of 
the ellipse in a fixed plane: variation of the ellipse amplitude 



n(t) in (a), the orientation angle 9(t) in (b), and the linearity 
X(t) in (c). These are modes of variability available to a 
modulated ellipse in two dimensions, as examined by |16|. 
The last two examples reflect new possibilities due to motion 
of the plane containing the ellipse: tilting of the plane due to 
variation of (3{t) in (d), and rotation of normal vector to the 
plane as a(t) varies in (e). This last mode can be visualized for 
a purely circular signal, X(t) = 0, as follows. Imagine a plate 
that is spinning on a table, with a particle running around 
the circumference of the plate. The spinning of the plate is 
associated with a(t), and as the plate slowly spins down f3(t) 
decreases to zero. 

As there are six parameters in the ellipse representation ©, 
and also six parameters in the Cartesian representation (0 
it appears reasonable to suppose that one set of parameters 
can be uniquely defined in terms of the other. While it is 
trivial to find expressions for the Cartesian amplitudes and 
phases in terms of the ellipse parameters, it is not so easy to 
accomplish the reverse. However, it is necessary to do so in 
order that the six ellipse parameters may be computed from the 
analytic versions of the three observed signals. This problem 
is addressed in the next two subsections. 

E. The Normal Vector 

A fundamental quantity, the normal vector to the plane 
containing the ellipse, will now be introduced. The normal 
vector n x (t) is defined as 

n x (t) = 3{x+(t)} x5R{x+(t)} (12) 

where "x" denotes the vector product or cross product and 
"3" the imaginary |>art0 That is, for two real-valued 3-vectors 
f = [fx fy fz] and g = [g x g y g z ] , "T" being the 
matrix transpose, the cross product is defined as 

f x g = Uy9z ~ fz9y)i 

- (fx9z - f z g x )j + (fx9y ~ f v 9x) k (13) 

where i, j, and k are the unit vectors along the x, y, and z- 
axes, respectively. Note that this definition of the normal vector 
n x (i) to the plane containing the motion is not the same as 
a more familiar quantity, the angular momentum vector; the 
relationship between these two quantities is outside the scope 
of the present paper and is left to a future work. 

With R a real orthogonal matrix such that R T = R _1 , and 
having unit determinant so that R is a proper rotation matrix, 
the cross product transforms as 

(Rf)x(Rg)=R(f xg) (14) 

a result that will be used repeatedly in what follows. This can 
be proven by first writing out the two sides in terms of the 
columns of R T , denoted [ri Y2 r^] = R T , which we note 
are related by e^k^k — *» X fj where e^-fc is the Levi-Civita 

4 The symbol "x" is also occasionally used herein to denote matrix 
multiplication or multiplication by a scalar, but the meaning will be clear 
from the context, since "x" can only denote a cross product when a vector 
multiplies another vector. 
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Fig. 2. Examples of modulated trivariate oscillations. In each one of the (a-e), only one of the five parameters describing the ellipse geometry varies. In 
(a-c), the plane containing the ellipse is held constant, but the ellipse amplitude n(t) varies in (a), the orientation 8(t) varies in (b), and the linearity X(t) 
varies in (c). The orientation of the plane containing the ellipse varies in (d) and (e), with the zenith angle /3(t) changing in (d) and the azimuth angle a(t) 
changing in (e). In each of these panels, a trace of the signal is shown in three dimensions with time visualized as the gray scale of the curve, with black 
being for early times. The projection of the signal onto the horizontal plane is also shown. The dashed line is a vertical line passing through the origin at 
x = 0, y = 0, while the plane z = is shown with a dotted line. The outline of the plane (or planes) containing the ellipse is also shown in each panel. Each 
of these signals is 800 points in length. The remainder of the caption pertains to panel (f), which is not referred to until Section [ill] In fact the five signals 
in (a-e) all have been constructed to have identical mean frequencies uJ x and mean second central moments <fj, as defined in Section [Hi] Panel (f) presents 
a spectral estimate of the joint spectrum Sx(ui) from each of the five trivariate signals, computed as described in the text. A dotted vertical line marks the 
global mean frequency uJ x /(27r) = 0.005 cycles per sample point or n X 10 — 2 radians per sample point. The spectra are virtually identical. 
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symbol. Equivalence between the two sides then follows from 
the Binet-Cauchy identity for four real-valued 3-vectors. 
The vector n x (t) may be expressed as, making use of ( TT4i >. 

n x (f)=J 3 (a(i))Ji(/3(t))J 3 (0(i)) 



a(t) sin (j)(t) 




a(t) cos cj)(t) 


-b(t) cos (j>(t) 


X 


b(t) sin </>(f ) 











= a(f)6(t)J 3 (a(f))Ji(/3(f))k (15) 

which is oriented perpendicular to the plane containing the 
ellipse and has magnitude ||n x (f)|| = a(t)b{t). Note 7r||n x (f)|| 
then gives the ellipse area. For future reference, we also define 



n x (t) = 1 j^W =J 3 ( a (t))J 1 ( / 9(t))k 



which is the unit normal to the plane containing the ellipse. 

Since the ellipse amplitude is already known through ©, 
there remain five ellipse parameters to solve for. From the nor- 
mal vector one may determine three further ellipse parameters. 
The linearity is found at once from 



A(f) 



a 2 (t) -b 2 (t) 



Mf)|| 



a»(t) + v H*+(*)H 4 ' 

Writing out components, the normal vector becomes 





sina(f) sin (3(t) 




n x {t) 


n x (i) = a(t)b(t) 


— cosa(f) sin /3(f) 




n y (t) 




cos /3(f) 




n z (t) 



(17) 



(18) 



and hence the angles (3(t) and a(t) may be readily determined. 
The former is 



/3(f) = 3 {in [n,(t) + iy/nl(t)+n*(t)\ } 
which recovers < /3(f) < tt, while the latter is 
a{t) = ${]n[-n y (t)+in x (t)]} 



(19) 



(20) 



giving — tt < a(f) < tt as desired. The use of the "9 {In [•]}" 
combination amounts to the so-called four quadrant inverse 
tangent function, with the usual choice that 3 {in [e ] } 
returns an angle 9 between — tt and tt. To see (T% , for example, 
note that 



5 {in \n z (t) + iy/n*(t)+nl(t)\ } = 

'{In [cos /3(f) + % |sin/3(f)|]} = 3? {lne^^} = /3(t) 



(21) 



substituting from ( TT8l on the second line, and using the fact 
that |sin/3(f)| = sin /3(f) since < /3(f) < tt by assumption. 

F. Precession and Phase Angles 

The values of four of the six ellipse parameters have now 
been established in terms of the canonical set of amplitudes 
and phases. To obtain the remaining two parameters, the 
orientation angle 6(t) and phase angle <j>(t), a representation is 
introduced that separates two-dimensional effects from three- 
dimensional effects in x + (f). 



Let H be the 3x2 matrix which projects a 2-vector onto 
the x—y plane in three dimensions, i.e. 



H 



Then one may write the analytic 3 -vector in (JSJ as 

x + (f) = [J 3 (a(f))J 1 (/3(f))H]x+(f) 



(22) 



(23) 



where x + (f) is a 2-vector which describes a modulated 
elliptical signal lying entirely within a plane. This complex- 
valued 2-vector is projected into three dimensions, tilted, and 
then rotated to generate the analytic 3-vector x + (f). Noting 



(16) [J 3 (a(f))J x (/3(f))H] T [J 3 (a(t))JiG8(t))H] 



H T H 



one can rearrange 

X+(f) : 



to find 



[J 3 (a(f))J 1 (/3(f))H] i x+(f). 



(24) 



(25) 



As a(t) and /3(f) have already been determined from the 
previous subsection, the 2-vector x + (f ) is now known for any 
given analytic 3-vector x + (f). 

The angles <p(t) and 9(t) may now be determined from 
x + (f), following |[T6l . Introducing the 2x2 counterclockwise 
rotation matrix 



cos ( 
sin( 



— sm( 
cos 9 



3(6) = 

the 2-vector x+(f) may be expressed as 

x+(i) = e^ (t) J(<9(f)) 



a(f) 
-ibit) 



(26) 



(27) 



which is the form for a modulated elliptical signal in two 
dimensions examined by lfl4l . ifTBI . Note that inserting (|27| | 
into (l23l gives ©, as required. Now define a new 2-vector 



.(f) 



a+(f)e^+W 


1 


'1 i ' 


o_(t)e^-W 


"71 


1 -i 



x+(f) (28) 



the amplitudes and phases of which are uniquely determined 
by the 2-vector x + (f). As discussed in IfTBI . z + (f) represents 
the motion in two dimensions in terms of the amplitudes and 
phases of counterclockwise-rotating and clockwise rotating 
circles, and leads to simpler expressions for the ellipse pa- 
rameters than does the use of x + (f). 

Substituting ( f27b for x+(f) into ( f28l . one finds z + (f) is 
expressed in terms of the ellipse parameters as 

- b{t)] e l9 ^ ' 



-(*) 



V2 



[a(t) 

[a(f) - 6(f)] e~ w ^ 



and so the orientation and phase angles of the ellipse are 



m = 

9(t) = 



4>+{t) + 4>-{t) 
4>+{t)-4>-{t) 



/2 
/2. 



(29) 

(30) 
(31) 



All six ellipse parameters are now uniquely determined in 
terms of a given analytic 3-vector x + (f). The functions a(f), 
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b(t), 6{t), 4>{t), a(t), and f3(t) so denned are called the 
canonical ellipse parameters. 

The 2-vector x + (i) describes the projection of elliptical mo- 
tion in three dimensions onto the plane which instantaneously 
contains the ellipse. A subtle point is that a(t), b(t), 6{t), 
and <fi(t) determined above are not necessarily the canonical 
ellipse parameters for this two-dimensional motion considered 
in isolation. This arises due to the fact that x + (t), and similarly 
z+(t), is not necessarily analytic. The message is that the 
canonical ellipse parameters give a unique description of the 
motion considered as a whole. This is analogous to the key 
point made by 132 '] for the univariate case that a x {t)e l ^ x ^ 
being analytic does not imply that e"^*) is also analytic. 

Note that choosing a different form for the representation of 
the modulated ellipse, using an alternate rotation convention 
such as x-y-z for example, would be equivalent to ©. The 
normal vector to the plane containing the ellipse, defined in 
(1121 . does not depend upon the particular ellipse representa- 
tion. Consequently in d23l one could have a different repre- 
sentation for the rotation matrix inside the square brackets, but 
its value must be unchanged. The parameters a{t), b(t), <f>{t), 
and 9(t), all of which are determined by the projection of the 
motion onto the plane instantaneously containing the ellipse, 
are thus also unchanged by an alternate rotation convention. 

III. Trivariate Instantaneous Moments 

Here the first- and second-order instantaneous moments of 
a trivariate signal are introduced and expressed in terms of the 
ellipse parameters. These time-varying quantities provide the 
link between the ellipse parameters and the Fourier spectrum 
of the signal. A fundamental quantity termed the trivariate 
instantaneous bandwidth is seen to capture five different modes 
of variability of ellipse geometry, all of which contribute to 
the second central moment of the signal's Fourier spectrum. 



A. Definitions 

This section will make use of the joint instantaneous 
moments of a multivariate signal introduced recently by [16]. 
These quantities integrate to the global moments of the aggre- 
gate spectrum of a multivariate signal, just as the instantaneous 
moments of a univariate signal integrate to the global moments 
of its spectrum [28|-[35|. The aggregate frequency-domain 
structure of the analytic vector x+(t) is described by the 
(deterministic) joint analytic spectrum 



5 x H^f- 1 ||x + H|| 



(32) 



where the total energy of the multivariate analytic signal is 
£ x =— / ||X + H|| 2 ^= / ||x+(t)|| 2 ^. (33) 

27r JO J-oo 

S x (uj) is the average of the spectra of the N analytic signals, 
normalized to unit energy. The joint global mean frequency 

1 f°° 

lj x =—J uS x (ui)du (34) 



is a measure of the average frequency content of the multi- 
variate analytic signal x + (i), while the joint global second 
central moment 

1 



cr x = — / (u - w x ) 2 5x(«) dw 
27T Jo 



(35) 



gives the spread of the average spectrum about the mean 
frequency. In the frequency-domain integrals above, the in- 
tegration begins at zero since X + (w) has vanishing support 
on negative frequencies by definition. 

The joint instantaneous frequency and joint instantaneous 
second central moment are then defined by lfl6l to be some 
time-varying quantities which decompose the corresponding 
global moments across time, i.e. which satisfy 



= £Z 



:|x+(t)|| 2 w x (t) dt 

\\x+(t)\\ 2 oi(t)dt 



(36) 



(37) 



noting that ||x + (i)|| 2 is aggregate instantaneous power of 
the analytic signal vector. Although the integrand in these 
expressions is non-unique, ifTBI show that the definitions 

3{x?(i)<(t)} 



UJ- 



:(*) 



!(*) 



I|x+(t)|| 
_(t) - v2 



-(*)! 



(38) 



(39) 



l|x + W|| 2 

are the natural generalizations of the standard univariate defini- 
tion of the instantaneous frequency [28] and the instantaneous 
second central moment [35 1. Note that d3~8l and d39l satisfy 
( |36l > and d37| i respectively. Also, d39l is nonnegative-definite, 
like the global moment to which it integrates. 

Thus uj x (t) and a x (t) can be said to give the instantaneous 
contributions to the mean Fourier frequency and second central 
moment, respectively, or equivalently, to partition the first two 
Fourier moments across time. The second-order instantaneous 
moment can alternately be expressed by defining the squared 
joint instantaneous bandwidth 



(40) 



which is that part of the instantaneous second central moment 
not accounted for by deviations of the instantaneous frequency 
from the global mean frequency. For a univariate signal 
x(t) = a x (t) cos 4> x (t), [16] shows that this definition gives 
Vx(t) = \a' x (t) / a x (t)\, the univariate instantaneous bandwidth 
identified by Il33l - ll35l . On account of the constraints d36l l and 
(|37| |. the scalar-valued functions ui x (t) and v x (t) summarize 
the time-varying frequency content of a multivariate signal, 
and its spread about this frequency, in the same manner in 
which the standard instantaneous frequency and bandwidth 
would accomplish this for a univariate signal. 

To find an expression for v x (t), we insert ( f4Qb into d39l to 
give, after some manipulation, 

II 1 1 2 

2 \\x' + (t) - iuj x (t)x+{t)\ l 



«£(*) 



(41) 



l|x + (*)ll 2 

which is the normalized departure of the rate of change of 
the vector-valued signal from a uniform complex rotation at a 
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single time-varying frequency w x (t). By contrast, <r x (f) i s by 
definition the normalized departure from a uniform complex 
rotation at a single fixed frequency ZU X . The squared multi- 
variate instantaneous bandwidth can alternately be expressed 

*+(*) 



(42) 



,|x + (f)|| 2 

in which we have made use of the definition of the multivariate 
instantaneous frequency ( l38l l. This form implies that when the 
modulus of the rate of change of x + (t) matches that expected 
for a set of sinusoids all locally progressing with frequency 
w x (f), the instantaneous bandwidth vanishes. 

For the trivariate case, it is desirable to obtain expression for 
the instantaneous moments in terms of the ellipse parameters. 
This would show how variations in the ellipse geometry 
contribute to the shape of the average spectrum, and at the 
same time provide a means for describing details of signal 
variation that are not resolved by the global moments. 

B. Trivariate Instantaneous Frequency and Bandwidth 

Forms for the trivariate instantaneous frequency and band- 
width in terms of the ellipse parameters will now be presented. 
The derivations of these forms are somewhat tedious and are 
therefore relegated to Appendix |B] here we will emphasize 
their interpretation. The trivariate instantaneous frequency is 

(i) (ii) 



w*(t) = w*(t) + y/i - A 2 (f) [w B (t) + u a (t) COS j8(t)] (43) 

and consists of two terms, (i) phase progression of the "par- 
ticle" along the ellipse periphery and (ii) the combined effect 
of internal precession and external precession of the ellipse. 
The squared trivariate bandwidth takes the form 



vl{t) 



(i) 
«'(*) 



«(*) 



(ii) 

41 - A 2 (f) 



(iii) 



(iv) 



+ A 2 (f ) [cue (t)+cu a (f ) cos f3(t)} 



(44) 



ll*+(*)H 

and consists of four terms, each of which is nonnegative: 
(i) amplitude modulation, (ii) deformation, (iii) precession, 
and (iv) the squared magnitude of that portion of the rate of 
change of the analytic signal that does not lie within the plane 
instantaneously containing the ellipse. 

If the motion is contained entirely within a single plane 
for all time, then a(t) and /3(f) are constants, and x' + (£) has 
no component parallel to the normal vector n x (f). Setting 
w Q (£) = and neglecting term (iv) in (l44l recovers the 
bivariate forms of the instantaneous frequency and bandwidth 

of ma 



(i) 



(ii) 



c(f) = LO lP (t) + y/l- \ 2 {t)uJ 9 (t) 



There are two new effects in the three-dimensional case 
compared with the two-dimensional case. Firstly, in both the 
trivariate instantaneous frequency and bandwidth, the effect of 
internal precession u>e(t) is modified by a term uj a (t) cos /3(f) 
which contains the external precession rate w Q (t). Secondly, in 
the trivariate instantaneous bandwidth (l44l . term (iv) emerges 
as a qualitatively new effect. 

Since the external precession uj a (t) is a frequency associ- 
ated with the azimuthal motion of the normal vector around 
the z-axis, it is clear that ui a (t) cos (3(t) is the component 
of external precession parallel to the normal of the plane 
containing the ellipse. Thus LJe(t) + u> a (t) cos f3(t) could be 
termed the "effective precession rate". When either (a) the 
plane containing the ellipse does not precess, so that uj a (t) 
vanishes, or else (b) the plane containing the ellipse rotates 
about a line it contains, so that (3(t) = tt/2, then the trivariate 
instantaneous frequency, and term (iii) in the bandwidth, both 
contain no contribution from the external precession. More 
generally, we may write the effective precession rate as 



LUg(t) + UJ a (t) cos /3(f) 



Vl-A 2 (i) 



(47) 



and in this form, term (iii) of the squared instantaneous 
bandwidth is identical in the bivariate and trivariate cases. 



C. Three-Dimensional Effects in the Bandwidth 

Term (iv) in d44l > involves the component of the rate of 
change of the analytic signal vector x + (f ) that lies parallel to 
the normal vector. That is, it is associated with the portion of 
changes in the real and imaginary parts of the analytic signal 
vector that deviate from the plane instantaneously containing 
the ellipse. In Appendix iBl this term is found to take the form 



n; 



-u; a (f) sin /3(f) wp(t)] x+(t) 



-(t)\\ 



(48) 



Thus changes in the zenith angle of the normal vector to the 
ellipse with respect to the vertical, i.e. the nutation rate U)p{t), 
as well as the component of precession of the plane containing 
the ellipse that lies in the horizontal plane, i.e. ui a (t) sin /3(f), 
both contribute. 

Writing out d48b in terms of the ellipse parameters leads to 
a rather complicated expression on account of the dependence 
of the 2-vector x+(£) on the orientation angle 6(t). However, 
one can apply the Cauchy-Schwarz inequality to d48b to find 



:(*K(t)i 

|x+(f)|| 2 



<w 2 (f)sin 2 ^(f)+a;?(f). 



(49) 



(45) 



An upper bound on the squared trivariate bandwidth is then 
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after applying the triangle inequality to term (iii) of the squared 
trivariate bandwidth, and setting A(t) to its maximum value 
of unity in this term. The bound (|50t has a considerably 
simpler form than (|44V The rates of change of each of the 
five parameters of the ellipse geometry appear. Note that the 
internal and external precession rates again contribute to the 
same nonnegative term. 

D. Invariance to Coordinate Rotations 

An important point is that the terms appearing in the 
expressions for the trivariate instantaneous frequency uj x (t) 
and squared instantaneous bandwidth v x (t) are independent of 
the choice of reference frame. That is, replacing the original 
real-valued trivariate vector x(t) with a rotated version, Rx(t) 
with dctR = 1, not only preserves the values of u x (t) and 
v x (t), as shown by [16], but also keeps the same values for the 
component terms. For example, the term w x (t) is determined 
from x + (t) and is independent of coordinate rotations, while 
w x (t), u)$(t), and u>$(t) are determined from x+(t) and thus 
are also invariant to rotations; then (l47t shows that the effective 
precession rate is likewise invariant. Hence all four terms in 
the squared trivariate bandwidth keep the same value under 
coordinate rotations. By contrast, the original definitions (|38l 
and (140b involve sums over terms in each signal component, 
and while the entire quantities u) x (t) and u x (£) are invariant to 
coordinate rotations, the contributing terms from each signal 
component are not. 

IV. Applications 

Two examples will be presented which illustrate the utility 
of the trivariate instantaneous moments. The first example 
returns to the synthetic signals constructed in Fig. |2] while 
the second examines a seismic record. 

A. Synthetic Example 

The five signals shown in Fig. |2^-e have been constructed 
such that the joint instantaneous frequency u x (t) has identical 
and constant values in each panel, as does the joint instanta- 
neous bandwidth |u x (i)|. Thus, the first two Fourier-domain 
moments of the aggregate spectra S x (ui) corresponding to 
these five signals should be identical apart from complications 
arising from the finite duration of the samples. The joint 
instantaneous frequency w x (t) takes a value of it x 10 -2 
radians per sample point, while |u x (t)| is set to 2.5tt x 10~ 4 
radians per sample point. The estimated aggregate spectra 
S x (lj) associated with each of the five trivariate signals are 
shown in Fig. [2J. These are formed with a standard multitaper 
approach [22|, [45 1 by averaging over three "eigenspectra" 
from data tapers having a time-bandwidth product P = 2; see 
P31 for details. 

It is seen that the five completely different rates of change 
all lead to nearly identical estimated spectra. The estimated 
aggregate spectrum therefore cannot distinguish between these 
different types of joint structure. However, the five rates of 
change can be directly computed from the observed signal, 
and therefore the different time-varying contributions to the 



Fourier bandwidth are known from the trivariate instantaneous 
moment analysis. A sixth possibility also exists. With the 
ellipse geometry fixed, we have ui x (t) = u>$(i) and therefore 
<r x (t) = [uj<j,(t) — uJ x ] 2 for the first instantaneous moment 
and second central instantaneous moment, respectively. Thus 
an appropriate uniformly increasing choice of oj^it) as the 
particle orbits a fixed ellipse could lead to a spectrum with 
the same first-order and second-order moments as in the five 
cases of time-varying ellipse geometry. 

As a caveat to this discussion, it is worth pointing out 
that for the short time series segments shown in Fig. |2^-e, 
the bandwidth of the estimated spectra in Fig. |2f are mostly 
due to the taper bandwidth and not to the signal bandwidth. 
Had we used much longer samples of these processes, we 
would have found that the five spectra differ in shape (and 
thus higher-order moments) despite having virtually identical 
first and second moments, by construction. 

B. Seismogram 

As a sample dataset, the seismic trace from the Feb. 9, 
1991, earthquake in the Solomon Islands as recorded at the 
Pasadena, California station (PAS), is presented in Fig. [3^. 
This dataset is useful as an illustration since the structure 
is quite simple, and since it has already been examined by 
other authors [23 1, 11271 . The data is available online from 
the Incorporated Research Institutions for Seismology (IRIS) 
using the WILBER data selection interfac^fl The source is 
located at 9.93° S, 159.14° E, while the station is located at 
34.15° N, 118.17° W. The bearing from the station to the 
source is 12.3° south of due west. The x, y, and z time 
series are rotated 12.3° clockwise about the vertical axis to 
form the radial-transverse-vertical records shown in Fig. [3^, 
with direction of the first (radial) axis pointing away from the 
seismic source. 

The distinct arrivals of two different types of surface waves 
are clearly visible in the time series: the Love wave is a 
linearly polarized oscillation in the transverse direction, while 
the Rayleigh wave is a roughly circularly polarized wave in the 
radial/vertical plane. Note that the Rayleigh wave is retrograde 
elliptical: particle paths in this wave move towards the source 
when they are vertically high and away from the source when 
they are vertically low. This is opposite from a gravity wave at 
a fluid interface, which undergoes prograde elliptical motion 
in the radial-vertical plane. 

Taking the analytic parts of these time series, the mul- 
tivariate instantaneous moments can be found at once. The 
ellipse amplitude n(t) and linearity X(t), as well as the joint 
instantaneous frequency w x (t), are shown in Fig. [3j>— c. The 
Love wave-dominated early portion of the record, between 
the two vertical lines, is clearly identified as being linearly 
polarized, while motion dominated by the Rayleigh wave 
after the second vertical line is associated with small linear- 
ity indicating slightly noncircular motion. The instantaneous 
frequency associated with both waves is observed to increase 
with time. A distinction between the frequency content of the 
two waves is also seen, with a sudden drop in instantaneous 

'http://www.iris.edu/dms/wilber.htm 
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Characteristics of a Seismic Record 




Fig. 3. A three-component seismogram from the Feb. 9, 1991, Solomon 
Islands earthquake as recorded in Pasadena, California. In (a), the seismic 
traces are shown in the radial-transverse-vertical coordinate system. Panels 
(b) and (c) show the ellipse amplitude re(t) and linearity A(t), respectively, as 
computed from the analytic signal. Finally (d) presents the joint instantaneous 
frequency u} x (t). The x-axis is time in seconds from the beginning of the 
record. Vertical lines mark the approximate locations of the arrival times of 
the Love wave and the Rayleigh wave. 

frequency after the Rayleigh wave arrival. Near the beginning 
and the end of the record, rapid fluctuations of the linearity 
and the instantaneous frequency are consequences of a low 
signal-to-noise ratio. 

A more informative presentation of the time-varying polar- 
ization is given in Fig.@] Here the coordinates are the standard 
Cartesian directions — East, North, and vertical. In Fig. |4&, the 
instantaneous orientation of the real-valued signal is visualized 



by plotting the values of unit vector x(t) = x(f)/||x(f)|| as 
points on the unit sphere. In Fig. the orientation of the unit 
normal vector n x (i) to the plane containing the signal and 
its Hilbert transform is similarly shown. During the Rayleigh 
wave, the unit normal vector offers a much more compact 
description of the signal. The direction of the signal vector 
x(t) oscillates throughout the radial/vertical plane, whereas 
the unit normal vector n x (i) is quite stable in the positive 
transverse direction. Note that there is not a comparable set of 
points in the negative transverse direction. This orientation of 
the unit normal vector indicates retrograde elliptical motion in 
the radial/vertical plane. Such an orientation is expected but 
is difficult to visualize from the raw time series. 

During the time period of the Love wave, the unit sig- 
nal vector x(f) oscillates between pointing in the positive 
transverse direction and the negative transverse direction. For 
such a one-dimensional signal, the plane containing the signal 
and its Hilbert transform is ill-defined; consequently, the 
direction of unit normal vector n x (i) is scattered, most likely 
meaninglessly, over the radial/vertical plane. This illustrates 
that the ellipse parameters should be interpreted carefully for 
signals that are nearly one-dimensional. 



C. Further Possibilities 

The trivariate instantaneous moments can be applied di- 
rectly to a time series to obtain useful information about 
the time-dependent polarization and evolution. This works as 
an initial analysis step when the signal is not too structured 
and when the signal-to-noise ratio is sufficiently large. These 
quantities can also be used as a building block in a more 
sophisticated analysis of realistic signals. It is now well known 
that examining the instantaneous moments of a composite 
univariate signal, consisting of the sum of several modulated 
oscillations, leads to unsatisfactory results [e.g. 06], and the 
same would be true for composite multivariate signals as 
well. In general, one would like to combine the trivariate 
instantaneous moments with a method for extracting different 
modulated oscillations from an observed signal vector. In the 
seismic example presented here, for example, one would like 
to form the trivariate instantaneous moments of the estimated 
Rayleigh wave signal and the estimated Love wave signal 
considered separately. 

The instantaneous moment analysis developed herein inter- 
faces well with the multivariate wavelet ridge analysis recently 
proposed by [47|. That method employs a time-frequency 
localization via the wavelet transform to reduce noise and to 
isolate different signal components from one another in a time- 
varying sense. The interface between these two methods is 
straightforward, as one can simply take an estimated analytic 
signal vector from the wavelet ridge analysis and determine 
its ellipse parameters from the model ©. As real-world data 
is generally noisy, in most applications of the instantaneous 
moment analysis it will be preferable to replace the true 
analytic signal with one estimated from wavelet ridge analysis 
or some related method. 



IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED AUGUST 201 1 



12 



Unit Signal Vector 



Unit Normal Vector 





-0.5 



1 -1 



] -1 



Y-Axis 



Y-Axis 



X-Axis 



X-Axis 



Fig. 4. Visualization of three-dimensional polarization of the Solomon Islands signal using (a) the real-valued unit signal vector x(t) = x(t)/||x(t)|| and (b) 
the unit normal vector n x (t). In both panels, the black "+" symbols show the vector position on the unit sphere during the Love wave event, i.e. between the 
two dotted lines in Fig. [3^, while the gray dots show the vector positions during the Rayleigh wave event, i.e. after the second dotted line. When appearing 
on the opposite side of the sphere, both types of symbols are plotted with a lighter shading. The usual Cartesian coordinate system is used. Heavy solid gray 
lines mark the positive x, y, and z coordinate axes, while heavy dashed gray lines mark the negative axes. A black circle marks the location of the negative 
transverse axis at 12.3° east of due south. 



V. Conclusions 

Any real-valued trivariate signal can be described as the 
trajectory traced out by a particle orbiting an ellipse, the 
amplitude, eccentricity, and three-dimensional orientation of 
which all evolve in time. The rate at which the particle 
orbits the ellipse, together with the rates of change of the 
ellipse geometry, control the first two moments of the Fourier 
spectrum of the signal. This perspective should be particularly 
valuable for describing signals which locally execute elliptical 
oscillations, but which may have broadband spectra on account 
of amplitude and frequency modulation — a class of signals 
which is expected to include many physical phenomena. 

While the link between instantaneous quantities derived 
from the analytic signal and the Fourier moments is well 
known for the standard univariate case, the instantaneous 
amplitude, frequency, and bandwidth take on geometrical in- 
terpretations in the trivariate case that enable a rich description 
of the signal's variability, permitting the distinction between 
qualitatively different types of motion. Compared with the 
bivariate case, new terms emerge in both the instantaneous 
frequency and bandwidth due to motions of the plane contain- 
ing the ellipse. As a consequence, there are six different ways 
that the spectrum of a trivariate signal, averaged over the three 
signal components, may have identical mean frequency and 
mean bandwidth, but arising from six very different pathways 
of evolution of the ellipse properties. 
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Appendix A 
A Freely Distributed Software Package 

All software associated with this paper is distributed 
as a part of a Matlab toolbox called Jlab, available at 
http://www.jmlilly.net The routines used here are mostly 
from the Jsignal and Jellipse modules of Jlab. The analytic 
part of a signal can be formed with anatrans. Given 
an analytic signal, instfreq constructs the instantaneous 
frequency and bandwidth, as well as the joint instantaneous 
frequency and bandwidth of a multivariate analytic signal. 
An elliptical signal in two or three dimensions is created by 
ellsig given the ellipse parameters, whereas ellparams 
recovers the ellipse parameters from a bivariate or trivariate 
analytic signal. The various ellipse geometry terms in the 
instantaneous bivariate or trivariate bandwidth are computed 
by ellband. Multitaper spectral analysis is implemented by 
mspec using data tapers calculated with sleptap. Ellipses 
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are plotted in two dimensions using ellipseplot. Finally, 
makef igs_trivariate generates all figures in this paper. 
All routines are well-commented and many have built-in 
automated tests or sample figures. 

Appendix B 

Expressions for the Instantaneous Moments 

In this appendix, the forms of the trivariate instantaneous 
frequency and bandwidth in terms of the ellipse parameters are 
derived. Some additional notation will facilitate the derivation. 
For a complex-valued 3-vector z(f), let 



as may be verified by direct computation. From ([TBI , we have 



*«(*) = [fi£(*M*)] fix(*) 

Z ± (t) = z(f)-Z||(t) 



(51) 
(52) 



be the components of z(f) instantaneously parallel to, and 
perpendicular to, the unit normal vector ii x (f) of the ellipse. 
The real and imaginary parts of zn(f) lie along the direction 
of n x (f), while the real and imaginary parts of z±(t) are in 
the plane perpendicular to n x (f). Note that the parallel part 
of the analytic signal vector vanishes, 

x|,(t)= [8£(t)R {*+(*)}] 8 x (f) 

+ i[n£(f)3{x+(f)}]n x (t) = (53) 

since by definition the unit normal is perpendicular to both the 
real and imaginary parts of x + (f); thus x+(f) = xj_(t). 

Using the parallel and perpendicular parts, we decompose 
the derivative of the analytic signal vector as 



[<(*)] || + [<(<)] X = X ||(*)+ X '-L( i ) 



(54) 



where we define the symbols xj,(f) and x' ± (t) to mean the 
parallel or perpendicular part of the derivative of x+(t). The 
action of taking the parallel or perpendicular part does not 
commute with the derivative, thus in general xj, (f), the parallel 
part of the derivative of x + (f), is not the same as the derivative 
of the parallel part of x + (f). The latter quantity is 

[x u (t)]' = xj,(t) + {[H x (f)] T x+(f)} n x (t) = (55) 

but since X|i(t) vanishes hence [xi|(i)]' does also, we have 

xI | (i) = -{[n x (t)] T x+(t)}n x (i) (56) 

as an expression for x'| (t) in terms of the rate of change of 
the unit normal vector n x (f). 

To find expressions for the rates of change of x+(f) and 
ii x (t) in terms of the ellipse parameters, note that 

[J 3 (a(i))J 1 (/3(t))]' = J 3 (a(i))Ji(0(t))x 





cos /3(t) 
-sin /3(f) 



-cos /3(f) sin /3(f) 



"0 ' 

-Wfi(t) 0-1 

1 



(57) 



n x (f)=J 3 (a(f))J 1 (^)) 



ui a (t) sin /3(f) 




(58) 



for the rate of change of the unit normal vector, which can 
have no component parallel to n x (f). Using (TS6b then gives 

xf,(t)= ([-w (t) sin /3(f) w /3 (f)] T x + (f))n x (f) (59) 

for the parallel component of the rate of change of x+(t). 
The perpendicular component of the rate of change of x + (f ) 
is found to be 

xl(f)=J 3 Kf))J 1 (^(f))Hx 

{x' + (f) +u a (t) cos /3(f) Jx+(f)} (60) 

where we let J = 3(tt/2) with no angle argument be the 
2x2 ninety degree rotation matrix; (l60l is obtained by 
differentiating x + (f) as expressed by d23l and then using 157) . 

To simplify the expression for x' ± (f), introduce for nota- 
tional convenience the two-vector 



r(f) 



" a(t) ' 


= «(*) 




A(f) 


-ib(t) 




- A(t)J 



(61) 



and then quadratic forms involving r(f) may be readily verified 



r H (f)r(f) - a 2 (f)+6 2 (f) 

r H (f)r*(f) = a 2 (f)-6 2 (f) 

r H (f)Jr(f) = 2ia(t)b(t) 

r H (f)Jr*(f) = 



which will be used shortly. Also we find 

A'(f) 



r'(f) = ^r(f) 



1 
I — 



Vl-A 2 (t) 



Jr*(t) 



(62) 
(63) 
(64) 
(65) 



(66) 



for the time derivative of r(f), using the definitions (0 and 
( flOb of the ellipse amplitude n(t) and linearity A(f). 

The rate of change of the two-vector x+(f) appearing in 
d60b is given by 



x' + (f) = e^ (t) J(0(f))x 



{r'(f) + i^(f)r(f) + w fl (t)Jr(t)} (67) 
as we find by differentiating (1271 1. Here we have made use of 
d 



dt 



J (0(f)) = J (6(t) + tt/2) = w fl (t) J (0(f)) J (68) 



for the derivative of the 2x2 rotation matrix. One finds 



i'(t) = e^«J(0(t)) x 



n{t) 



iuj<t,(t) 



r(f) 

Jr*(f)l (69) 



after making use of d66l i for r'(f). 
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In terms of the parallel and perpendicular components of [10 



x^_(i), the instantaneous frequency and bandwidth become 



*(*) = 



;(*) = 



«(tK(f)} 



l|x + (*)|| 

«ll 2 



-(*)« 



-(*)! 



(70) 



-<*£(*) (71) 



using (l40l for the latter, and noting x^(i)xj,(t) = 0. For the 
instantaneous frequency, substituting ( f60b into (T70b gives 

x?(t)xl(t) 

= x£ (t)x' + (t) + Lo a (t) cos f3(t)*+ (t) Jx+ (t) (72) 

and using (|62l-(|o3]> together with ( l67l i, the trivariate instanta- 
neous frequency expression ( |43l then follows. For the trivariate 
instantaneous bandwidth, note that the first term on the right- 
hand-side of (Tm > becomes, substituting (f60b for Xj_(t), 



.(*)« 



.(*)! 



|x+(t)|| 



We find using 



llx+WII 2 
+ 2w a (t) cos 13 (t) 

that 



uj 2 a {t) cos 2 f3(t) 

R{x?(t)J T xV(t)} 



-(*)ll 



din 



1-A 2 (t) 



■(*)« 



1 dA(t) 
2 



+ (t) + 2 Vl- A 2 (<) w^ (t) we (t) + w| (t) (74) 



and similarly 

SR{xf(*)J r x' + ft)} 
l|x + (i)|| 2 



w e (t) + v/l-A 2 (t)w^(i). (75) 



Combining ( 174b and ( fTBI l into d73"l l. then using this together 
with d43l for w x (i) in ( PTTl i. cancelations occur, leading to the 
form of the trivariate bandwidth (l44l given in the text. 
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